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Abstract 

In the first part of this paper, we propose new optimization-based 
methods for the computation of preferred (dense, sparse, reversible, 
detailed and complex balanced) linearly conjugate reaction network 
structures with mass action dynamics. The developed methods are ex- 
tensions of previously published results on dynamically equivalent reac- 
tion networks and are based on mixed-integer linear programming. As 
related theoretical contributions we show that (i) dense linearly conju- 
gate networks define a unique super-structure for any positive diagonal 
state transformation if the set of chemical complexes is given, and (ii) 
the existence of linearly conjugate detailed balanced and complex bal- 
anced networks do not depend on the selection of equilibrium points. 
In the second part of the paper it is shown that determining dynami- 
cally equivalent realizations to a network that is structurally fixed but 
parametrically not can also be written and solved as a mixed-integer 
linear programming problem. Several examples illustrate the presented 
computation methods. 

Keywords: chemical kinetics; stability theory; weak reversibility; linear 

programming; dynamical equivalence 

AMS Subject Classifications: 80A30, 90C35. 



1 Introduction 



The mathematical study of chemical reaction networks is a rapidly grow- 
ing field that has been applied recently to research problems in industrial 
chemistry, systems biology, gene regulation, and general nonlinear systems 
theory, among others 10 23|[33] . There has also been significant theoretical 
work in the literature on such topics as persistence [T}|3| [26|[32] , multistabil- 
ity [8|[9|[3l], monotonicity [5|[6], the global attractor conjecture for complex 
balanced systems Jl| 2|[7|[! 



tion networks 11,25 



lumping [15 ,28 39 40^, and conjugacy of reac- 
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One line of research which has been attracting increased attention has 
been that of determining when two reaction networks exhibit the same qual- 
itative dynamics despite disparate network structure. In [ll] and 34 , the 
authors complete the question of what network structures can give rise to 
the same set of governing differential equations and therefore exhibit iden- 



tical dynamics. This work was extended in 25 to networks which do not 
necessarily have the same set of differential equations but rather have tra- 
jectories related by a non-trivial linear transformation. Similar ground has 
also been touched in the papers 28 and [Ts] which deal with properties of 
linear lumpings, linear mappings which potentially reduce the dimension of 
the species set. 

In this paper we look at the problem of algorithmically determining 
when a network is linearly conjugate to another network satisfying speci- 
fied conditions. This problem was first addressed in [35] where the author 
presents a mixed-integer linear programming (MILP) algorithm capable of 
determining sparse and dense realizations, i.e. networks with the fewest and 
greatest number of reactions capable of generating the given dynamics. The 
algorithm was extended to complex and detailed balanced networks in [36| , 
reversible networks in [37| , weakly reversible networks in [38] , and linear con- 
jugate networks in [27], where a computationally more efficient procedure 
for determining weak reversibility was also presented. 

We will continue in this paper the development and application of this 
MILP framework to linearly conjugate networks. In particular, we will give 
conditions which can be used to find sparse and dense linearly conjugate 
networks which are detailed and complex balanced, fully reversible, and 
which contain the greatest and fewest number of complexes. We will also 
expand the original MILP algorithm for finding sparse and dense realizations 
to find alternative realizations to a given reaction network when the network 
structure is fixed but the parameter values are not. Since this algorithm 
works without having to have the rate constants specified beforehand, this 
will allow us to answer questions about the reaction mechanism itself. 



2 Background 

In this section we present the terminology and notation relevant to chemical 
reaction networks and the main results from the literature upon which we 
will be building. 
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2.1 Chemical Reaction Networks 

We will consider the chemical species or reactants of a network to be given 
by the set S = {Xi,X2, ■ ■ ■ , Xn}. The combined elements on the left-hand 
and right-hand side of a reaction are given by linear combinations of these 
species. These combined terms are called complexes and will be denoted by 
the set C = {Ci, C2, . . . , Cm} where 

n 

Ci = '^aijXj, i = l,...,m 
i=i 

and the Oij are nonnegative integers called the stoichiometric coefficients. 
We define the reaction set to be 7^ = {(Q, Cj) \ Ci reacts to form Cj} where 
the property {Ci,Cj) G TZ will more commonly be denoted Ci — )■ Cj. To 
each {Ci,Cj) £ TZ we will associate a positive rate constant k{i,j) > and 
to each (Ci,Cj) 7^ we will set k{i,j) = 0. The triplet J\f = {S,C,n) wih 
be called the chemical reaction network. 

The above formulation naturally gives rise to a directed graph C{V, E) 
where the set of vertices is given hy V = the set of directed edges is given 
hy E = IZ, and the rate constants k{i,j) correspond to the weights of the 
edges from Ci to Cj. In the literature this has been termed the reaction 
graph of the network |24| . Since complexes may be involved in more than 
one reaction, as a product or a reactant, there is further graph theory we 
may consider. A linkage class is a maximally connected set of complexes, 
that is to say, two complexes are in the same linkage class if and only if there 
is a sequence of reactions in the reaction graph (of either direction) which 
connects them. A reaction network is called reversible if Ci — t- Cj for any 
Ci,Cj £ C implies Cj Ci. A reaction network is called weakly reversible 
if Ci — )• Cj for any Ci,Cj G C implies there is some sequence of complexes 
such that Ci = C^(i) C^(2) > C'^(j-i) C'^(i) = Cj. 

A directed graph is called strongly connected if there exists a directed 
path from each vertex to every other vertex. A strongly connected component 
of a directed graph is a maximal set of vertices for which paths exists from 
each vertex in the set to every other vertex in the set. For a weakly reversible 
network, the linkage classes clearly correspond to the strongly connected 
components of the reaction graph. 

Assuming mass-action kinetics, the dynamics of the specie concentra- 
tions over time is governed by the set of differential equations 

- = y.^,.f(x) (1) 
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where x = [xi 3:2 • • • x„]^ is the vector of reactant concentrations. The 
stoichiometric matrix Y contains entries [Y]ij = aji and the Kirchhoff or 
kinetics matrix is given by 

When we speak of the structure of a kinetics matrix, we wih be referring 
to the distribution of positive and zero entries, which determines the net- 
work structure of the corresponding reaction graph. FinaUy, the mass-action 
vector ^'(x) is given by 

n 

vi/,(x)=n^r'", j=i,...,m. (3) 



2.2 Linearly Conjugate Networks 

Under the assumption of mass-action kinetics, it is possible for the trajecto- 
ries of two reaction networks M and M' to be related by a linear transforma- 
tion and therefore share many of the same qualitative features (e.g. num- 
ber and stability of equilibria, persistence/extinction of species, dimensions 
of invariant spaces, etc.). This phenomenon was termed linear conjugacy 
in [25]. 

For completeness, we include the formal definition of linear conjugacy as 
presented in [25]. We will let <I>(xo,t) denote the flow of ([T]) associated with 
Af and ^'(xo,t) denote the flow of ([T]) associated with A/"'. 

Definition 2.1. We will say two chemical reaction networks N andN' are 
linearly conjugate if there exists a linear function h : M"q i— q such 
that h($(xo,t)) = ^'(h(xo),t) for all xq G M>o- 

It is known that linear transformations h : M" q i— t- M" q can consist of at most 
positive scaling and reindexing of coordinates (Lemma 3.1, |25j). Linear 
conjugacy has been subsequently studied from a computational point of 
view in [27| . 

Linear conjugacy is a generalization of the concept of dynamical equiv- 
alence whereby two reaction networks with different topological network 
structure can generate the same exact set of differential equations ([T]). Two 
dynamically equivalent networks N and M' are said to be alternate realiza- 
tions of the kinetics ([T]) , although it is sometimes preferable to say that N is 
an alternative realization of N' or vice- versa. Since the case of two networks 
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being realizations of the same kinetics is encompassed as a special case of 
linear conjugacy taking the transformation to be the identity, we will focus 
on linearly conjugate networks. 

In general practice, we are given a network M and asked to determine 
its dynamical behaviour. This is often a challenging problem; however, 
we may notice that the network M behaves like one from a well-studied 
class of networks and therefore suspect a relationship which preserves key 
qualitative aspects of the dynamics. The theory of linear conjugacy can 
provide a powerful tool in analyzing such networks. If the network can be 
shown to be linearly conjugate to a network M' from the class of networks 
with understood dynamics, the dynamics of M' are transferred to A^. 

This raises the question of how to find a linearly conjugate network J\f' 



when only the original network J\f is given. This was studied in 27 where 



the authors built upon the linear programming algorithm introduced in [35] . 
We can impose that a network M' be linearly conjugate to our given network 
M with the set of linear constraints 



(LC) { 



Y ■Ab = T-^ -M 

m 

^[Ab\ij=0, j = l,...,m 

[Ablij > 0, i,j = l,...,m,i^j 
[Ablii < 0, i = l,... ,m 
e<Cj<l/e, j = l,...,n 



where < e ^ 1, and the matrices M £ 



and T G 



(4) 



are given by: 



M 
T 



Y ■ Ai;., and 
diag{c} . 



(5) 
(6) 



The kinetics matrix for the network J\f' can by constructed from Aj, £ 



imxm J, g ^j^g relation 



A'k = Ab-diag{^{c)} . 



(7) 



Finding a network satisfying Q and then solving ([T]) is sufficient to deter- 
mine a network M' which is linearly conjugate to J\f via the transformation 
h(x) = T-^x. 



2.3 Sparse and Dense Linearly Conjugate Networks 

In order to place the problem within a linear programming framework, we 
need to choose an objective function to optimize. An appropriate choice 
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of such a function is not obvious and may vary depending on the intended 
apphcation. 

One particularly intuitive choice, which was introduced in ]35] and has 
been widely used since, is to search for networks J\f' with the fewest and 
greatest number of reactions {sparse and dense networks, respectively). A 
sparse (respectively, dense) linearly conjugate network is given by a ma- 
trix A'f^ satisfying Q with the most (respectively, least) off-diagonal entries 
which are zeroes. Since the structure of A'f^ and Af, are the same, a corre- 
spondence between the non-zero off-diagonal entries in A'j^ and a positive 
integer value can be made by considering the binary variables 5ij G {0, 1} 
which will keep track of whether a reaction is 'on' or 'off', i.e. we have 

6ij = I [Ab\ij > e, i,j = l,...,m, i^j 

where < e ^ 1 is sufficiently small and can be chosen the same as in 
Q, and where the symbol '<^=^' denotes the logical relation 'if and only 
if. These proposition logic constraints for the structure of a network can 
then be formulated as the following linear mixed-integer constraints (see, for 
example, [30j): 

<[Ab]ij - edij, i,j = l,...,m, i j 
< -[Ab]ij +Uij6ij, i,j = l,...,m, i^j (8) 
Sije {0,1}, hj = !,■■■, m,i^j, 

where Uij > for i,j = l,...,m,i ^ j, are appropriate upper bounds 
for the reaction rate coefficients. The number of reactions present in the 
network corresponding to ^4^ is then given by the sum of the dij^s so that 
the problem of determining a sparse network corresponds to satisfying the 
objective function 



(Sparse) 



minimize 



over the constraint sets Q and ([s]) . Finding a dense network corresponds to 
maximizing the same function, which can also be stated as a minimization 
problem as 

{m 
minimize (10) 

It is known that, for trivial linear conjugacies, the structure of the dense 
realization contains the structure of all other trivial linear conjugacies as a 
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subset (Theorem 3.1, [3^)- We now prove the comparable result for non- 
trivial linear conjugacies. 

Theorem 2.1. Let M he a chemical reaction network. Suppose that the 
reaction network M' is linearly conjugate to N and dense. Suppose that J\f" 
is also linearly conjugate to N . Then the directed unweighted graph of N" 
is a subset of the directed unweighted graph of N' . 

Proof. Assume N' and N" are linearly conjugate to N , N' is dense in the 
space of networks which are linearly conjugate to M , and M" contains a 
reaction not contained in M' . 

Since both M' and M" are linearly conjugate to A^, we have by Q that 

Y ■Ak = T' -Y ■A'i, 

and 

Y ■ Ak = T" - Y ■ A'l 

where T = diag{c'}, c' G M"q and A'^ correspond to TV', c" G M"o and 
T" = diag{c"} and A'l correspond to M" . 
Now consider T' - Y ■ A'l. We have 

T' -Y ■ A'l = T' ■ {T")-^ -T" ■Y-A'l = Q-Y-Ak 

where Q = T' ■ {T")^^. Consequently, we have 

T' -Y ■ A'l + T' - Y ■ A't^ = {Q + I) -Y ■ Ak. (11) 

On the other hand, we have 

T' -Y ■ A'l + T' - Y ■ A't, = T' -Y ■ « + A'f^) =T' -Y ■ A!l' (12) 



where A'l' = A'^^ + A'l. Combining (11) and (12) gives 

Y ■Ak = {Q + -T' -Y ■ A'l' = T'" ■ Y ■ A'l' 

where T'" = {T' ■ {T")-^ + I)-^-T'. 

Since T', T" and / are diagonal matrices with positive entries on the 
diagonal, so is T'" . This means that the network M'" corresponding to 
A'l' = A'^ + A'l is linearly conjugate to M . We can readily see, however, 
that M'" contains all of the reactions in both M' and N" . If N" contains a 
reaction not contained in N' then N'" clearly has more reactions than M' 
which contradicts the assumption that M' is dense in the space of networks 
which are linearly conjugate to N . The result follows. □ 
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The following result follows immediately. 



Corollary 2.1. Let M be a chemical reaction network. Then the structure 
of the unweighted directly graph of the dense reaction network M' which is 
linearly conjugate to M is unique. 



Proof. This follows directly from Theorem 2.1 



□ 



3 Computational Extensions of Linearly Conju- 
gate Networks 

In this section we will extend the optimization framework introduced in 



Section 2^ to include complex balanced, reversible and detailed balanced 
networks, and to search for networks with the greatest and fewest number 
of complexes. 



3.1 Weakly Reversible Networks 

Weakly reversible networks are a particular important class of reaction net- 
works because strong properties are known about their dynamics and equi- 
librium concentrations. Under a supplemental condition on the rate con- 
stants, weakly reversible networks are known to have complex balanced 
equilibrium concentrations and therefore exhibit all of the dynamical proper- 



ties normally reserved for these networks [17|22| (see Section 3.3 for further 
discussion of complex balanced networks). 

Consequently, they are a primary candidate for the type of network we 
would like to find. The problem of determining if and when a chemical 
reaction network J\f is linearly conjugate to a weakly reversible network M' 
was first considered in 38 and further refined in 27 . For convenience, 



we briefiy recall the constraints published in [27] that guarantee the weak 
reversibility of the linearly conjugate network N': 



(WR) < 



[Ak\ij > 0, i,j = l,...,m, iy^j 



1, . . . ,m 



(13) 



where is an auxiliary Kirchhoff matrix with the same structure as with 
appropriately scaled columns such that its kernel contains the m-dimensional 
vector of all ones. In order to guarantee that the matrix has the same 
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structure as A'j^ and we also require that 



(WR-S) 



l,...,m, ij^j 



Q < -[Ak]ij + UijSij, i,j = l,...,m, i^j. 



3.2 Reversible Networks 



In [36] , an algorithm was presented which was capable of determining re- 
versible reaction networks which are trivially linearly conjugate to a given 
reaction network. In this section, we present a simplified methodology and 
apply it to non-trivial linear conjugacies. 

We recall that a network is reversible if Ci — )• Cj for any Cj, Cj € C implies 
Cj —7- d. For the network M' , this is equivalent to the condition 

[A'khj > e ^ [A'k]ji > e 
for some sufficient small < e ^ 1. This is in turn equivalent to 



1 



1 



where 5ij £ {0, 1}, i,j = 1, . . . ,m, i ^ j, as in Section 2.3 It follows that 
we can restrict our search space to reversible networks with the constraint 
set 



(Rev) 



(15) 



1, . . . ,m, i < j. 

A sparse or dense reversible network which is linearly conjugate to J\f can 
be found by optimizing ^ or (10), respectively, over the constraint sets Q, 
P and 



3.3 Complex Balanced Systems 

A particularly important class of chemical reaction networks are the complex 
balanced networks introduced in [24| . 

Definition 3.1. An equilibrium concentration x* € M"o '^f chemical 
reaction network M is a complex balanced equilibrium concentration 
if 

Ak ■ *(x*) = 0. (16) 

The network Af is called complex balanced if every equilibrium concentra- 
tion X* G M!^Q is a complex balanced equilibrium concentration. 
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Many strong properties are known about complex balanced networks. In 
particular, it is known that complex balanced networks permit exactly one 
positive equilibrium concentration in each invariant space of the network and 
that this equilibrium concentration is locally asymptotically stable (Lemma 
4C and Theorem 6A, [24|). Complex balanced systems are also known to be 
weakly reversible so that they are a subset of the weakly reversible networks 



considered in Section 3.1 (Theorem 2B, [22]). 

The following result shows complex balancing is a system property de- 
pending on the structure and parameters of the network and not on the 
chosen equilibrium concentration. 

Theorem 3.1 (Theorem 6A, [24] ). If a chemical reaction network Af is 
complex balanced at an equilibrium concentration x* E M" q then it is complex 
balanced at all of its equilibrium concentrations. 

It should be noted that the complex balancing of a network is still dependent 
on the choice of rate constants. It is possible for a reaction network to be 
complex balanced for some choices of rate constants and not for others. 

In [36], an algorithm was presented which was capable of determining 
sparse and dense complex balanced networks which are trivially linearly 
conjugate to a given network J\f. This method required determining an 
equilibrium value x* G M" q of the network M and then imposing the condi- 
tion 

A', . M/(x*) = 



on M' in accordance with (16). In this section, we extend these results to 
include non-trivial linearly conjugate networks. 

Suppose that M and J\f' are linearly conjugate via the transformation 
y = T^^x. In order to guarantee the network A/"' is complex balanced. 



according to (16) we require that 

A', ■ *(y*) = 0. (17) 
Since the equilibrium concentrations of N and M' are related by the transfor- 



mation y* = T ^x*, we have that the left-hand side of (17) can be rewritten 

Ai . *(y*) = Ai . M/(r-ix*) = A', . diag{^{c)}-' ■ ^(x*) = A, • M/(x*) 

where we have made use of the form of the kinetics matrix of M' according to 
([7|. The condition for complex balancing of the linearly conjugate network 
M' is therefore 

Ab ■ ^'(x*) = 
(CB) { M-^'(x*) = (18) 



>0 
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where M is as in ([s]). A sparse or dense complex balanced network which is 
linearly conjugate to Af can be found by optimizing ([9|) or (10), respectively, 
over the constraint sets Q, ([8| and (18). 



It should be noted that the optimization algorithm is less computation- 



ally exhausting than the corresponding algorithm for weak reversibility ( 13 ) . 
This is because the matrix Af^ required in the general weakly reversible case 
is not required in the complex balancing condition; rather, it is sufficient 
to use the matrix Af,. Consequently, there are fewer decision variables in 
the complex balancing algorithm. The pre-step that a x* G M"o be found 
satisfying Y ■ A^ - ^'(x*) = may off-set this advantage, however, depending 
on the difficulty in solving Y ■ Aj. ■ ^'(x*) = 0. 

It is unclear how the outcome of the algorithm depends on the choice of 
equilibrium concentration x* E IK>0) which in general is not unique. The 
following result clarifies this dependence. 

Theorem 3.2. Suppose N is linearly conjugate to N' with transformation 
matrix T = diag{c} where c £ M"q and suppose N' is complex balanced at 
y* = r~^x* where x* E M"o '^^^ Y ■ Ak ■ '^{x*) = 0. Then N' is complex 
balanced at y* = T~^x* for all x* E M"o satisfying Y ■ A^ ■ ^'(x*) = 0. 

Proof. Suppose trajectories of M' are related to trajectories of M by the 
relationship y = T~^x where T = diagjc} and c E M"q. 

Suppose that M' is complex balanced at y* = T^^x* where x* E M?q is 



an equilibrium concentration of M. It follows from Theorem 3.1 that 



Y.A',.^{y) = =^ A',-^iy) = 0. (19) 

Now consider an arbitrary equilibrium concentration x* E M"q of M. 
Since J\f and M' are linearly conjugate, it follows that Y ■ Ak = T ■ Y ■ A'/^ ■ 
diag {^(c)}~^. It follows that we have 

= Y-Ak- ^{5t*) =T-Y-A'f,- diag{^(c)}-^ • ^(Ty*) 
= T-Y-A',-^{r). 



It follows by the structure of T that Y ■ A'^- ^'(y*) = 0. From (19) we have 
that A'j^ ■ ^(y*) = 0. In other words, TV' is complex balanced at y* and we 
are done. 

□ 

This result shows that when imposing the complex balancing constraint 



(18) on M' , it does not matter which equilibrium concentration of N we 



choose. The feasible set of solutions (i.e. admissible networks) is the same. 
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3.4 Detailed Balanced Systems 



In 36 , the authors present an algorithm for determining detailed balanced 
networks which are trivially linearly conjugate to a given network. In this 
section we extend this algorithm to non-trivial linear conjugacies. 

Definition 3.2. An equilibrium concentration x* G M"q of the chemical 
reaction network N is a detailed balanced equilibrium concentration if 

[Ak]ij-^j{-K*) = [Ak]ji-^,{^*), Vi,j = l,...,m, i^j. (20) 

The network J\f is called detailed balanced if every equilibrium concentra- 
tion X* G M"q is a complex balanced equilibrium concentration. 

In other words, an equilibrium concentration x* E M"q is detailed balanced 
if the flow across each reaction is balanced by the flow across an opposite 
reaction at x*. 

Suppose that Af and M' are linearly conjugate via the transformation 
y = T~^x. In order to guarantee the network J\f' is detailed balanced, 
according to (20) we require that 

diag{M/(y*)} • {A',f = A', • diag {vl/(y*)} . (21) 

Since the equilibrium concentrations of M and M' are related by the trans- 
formation y* = T~^x*, we have that 

A', ■ diag{^(y*)} = A', ■ diag {vl/(c)}-^ • diag{^(x*)} 
= A-diag{*(x*)} 

where we have made use of the form of the kinetics matrix of M' according to 
([T]). The condition for detailed balancing of the linearly conjugate network 
M' is therefore 

r diag{vI/(x*)}-Af = A-diag{M'(x*)} 
(DB) I M ■ ^'(x*) = (22) 

where M is as in ([5]). A sparse or dense detailed balanced network which is 
linearly conjugate to M can be found by optimizing ^ or (10), respectively, 
over the constraint sets (Q, ([s]) and (22). Since the analogous result to 
Theorem |3.2| holds as a consequence of detailed balanced systems being a 
subset of complex balanced networks, we do not prove it here. 
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3.5 Minimal and Maximal Number of Complexes 

We can also adapt to non-trivial linear conjugacies the algorithm introduced 
in [36] for determining a network with the fewest or greatest number of 
complexes from a fixed complex set which is trivially linearly conjugate to 
a given network A^. 

In order to count the number of complexes in the network, we introduce 
the binary variables 5i G {0,1}, i,j = l,...,m, and consider the logical 
equations 



5i = l ^ J2 i^khn + Yl [^-tlj^^ > 



(23) 



ji=i 



J2 = l 



for i = 1, . . . ,m. In other words, 6i takes on the value of one if and only if 
there is a reaction to or from the complex Ci in the network; otherwise, it 
takes the value zero. For computational purposes, we reconsider (23) as 



(24) 



ii=i 



J2 = l 



where < e ^ 1. The linear constraints required to substantiate (24) are 



(Comp) < 



ii=i 



i2=i 



0< - J2^Mm - J2^Ak]j2i 



ii=i 



( 



i2=i 
i2^« 



+ 



\ 



V2i 



.il = l 



i2=i 



5i e {0,1} ,i = 1, ... ,m. 



(25) 



We can now determine a network with the fewest or greatest number of 
complexes by optimizing the functions 



(Min) 



mmmiize 



i=l 



(26) 
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or 

{m 
minimize — 6i (27) 
i=i 

respectively over the constraint sets Q and (25). Further constraints can 
be imposed to restrict ourselves to specific classes of systems (e.g. complex 
balanced systems, reversible networks, etc.) although care has to be taken 
to ensure the structural constraints are still satisfied. 

3.6 Examples 

In this section, we present a few examples which illustrate the methodologies 
outlined so far. 

Example 1: Consider the kinetic scheme 

Xl = X1X2 — 2x\ + xix^ 

X2 = —x\x\ + X1X3 (28) 

2 Q 2 
x^ = Xl — 0X1X3 

introduced in j27j. In that paper, it was shown that the kinetics could be 
generated by a reaction network involving the complex set 

Ci = Xi + 2X2, C2 = 2X1 + 2X2, = 2X1 + X2, 
C4 = 2X1, = Xl, Ce = 2X1 + X3, C7 = Xl + 2X3 
Cs = 2X1 + 2X3, Cg = Xl + X2 + 2X3, Cio = Xl + X3 



and that (28) has dynamics which is linearly conjugate to those generated 
by the sparse weakly reversible network given in Figure [l][a) (conjugacy 
constants ci = 20, C2 = 2, C3 = 5) and the dense weakly reversible network 
given in Figure [l]||b) (conjugacy constants ci = 20/3, C2 = 20/33, C3 = 5/3). 

The network in Figure [l]^a) is complex balanced as a consequence of it 
being a zero deficiency weakly reversible network. It can be easily checked, 
however, that the network in Figure[T]^b) is not complex balanced. We might 
wonder, therefore, what running the algorithm for a dense complex balanced 



network which is linearly conjugate to a network generating the kinetics (28) 
would produce. 

Numerically, we can determine that an equilibrium concentration of 
([28j) is {x\,xl,x%) = (0.2,0.577350269,0.258198889). Running GLPK for 
a sparse network ^ over the constraints (Q , ([s]) and ([is]) gives the network 
given in Figure [T][c) (conjugacy constants ci = 20/3, C2 = 20/33, C3 = 5/3). 
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(a) X1+2X2 - 

125\^ 
X1+2X3 



25 



40 



2X1+2X2 

2x1 



(b) X1+2X2 2X^+2X, 



^2 

0.926^ 



13.9 



X1+2X3 <=^ 2X1 
1 ^ ^13.3 ' 

0.926 




2X1 +X2 



(c) X.+2X, 5^ 2Xi+2X 



13.9 



^2 
0.926^ 



12.5 



1.35 



2 
,2.01 



X1+2X3 



-U926: 



■13.3 



2>C|+X2 



0.926 



Figure 1: Weakly reversible networks which are linearly conjugate to a net- 



work with the kinetics (28). The network in (a) is sparse while the networks 



in (b) and (c) are dense. The networks (a) and (c) are also complex bal- 
anced. (The parameter values in (b) and (c) have been rounded to three 
significant figures.) 



We notice that this network has the same structure as the weakly reversible 
network in Figure [l]l|b) and, furthermore, only differs in two rate constant 
values. It can be checked, however, that (c) is complex balanced while (b) 
is not. 



Example 2: Consider the kinetic scheme 



xi = —2xiX2 + 2x3 + 2X6 
±2 = -X1X2 + 2x3 

X3 = 2xiX2 - 4X3 



X4 

Xq 



X3 - X4X5 + X6 

—2x4X5 + 4X6 

X4X5 — 2X6. 



(29) 



Using the indexing scheme introduced in 20 and more recently applied 



in the papers [36| and 27 we can construct a chemical reaction network 
capable of generating the dynamics (29) under the assumption of mass- 
action kinetics ([T]) which involves the complexes 

Ci = Xi+ X2, C2 = X2, C3 = X3, C4 = Xi + X3, 
C5 = Xq, Cq = Xi + Xq, C7 = Xi,Cs = X2 + Xs, 
Cg = Xi + X2 + X^, Clo = 0, Cii = X3 + X4 
C12 = X4 + X^, Ci3 = X5, Ci4 = Xq, Ci5 = X4 + Xq, 
CiQ = X4, Cn = X5 + Xq, Cis = X4 + X5 + Xq. 
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It seems less than desirable, however, to consider a network of 18 com- 



plexes given the simplistic dynamics (29). We might wonder if there is a 
network with fewer complexes. Optimizing (26) over the constraint sets (|4]) 
and (25) gives the network 

1 1 2 

X2 + ^3 < Xl + X2 ^ X3 7- X4 



2 ^2 
X5 i — X4 + ^ Xq — > Xl + 

2 



and the conjugacy constants ci = 1, C2 = 2, C3 = 2, C4 = 1, C5 = 4, ce = 2. 



In terms of understanding the qualitative dynamics of 29 ) , this network 
is not particularly insightful. We might notice, however, that the net effect 
of all the reaction pathways leading out from X3 is to create an X2 and an X4 
at the expense of depleting X^. The complex X2+X4, however, has not been 
considered in the procedure. Similarly, the reaction pathways leading out 
from Xq generate an Xi and an X^ but the complex Xi + X^ has not been 
included in the procedure. We might consider appending the procedure, 
therefore, to include C19 = X2 + X4 and C20 = Xi + X^. Repeating the 
algorithm in GLPK gives the network 



1 



2 



Xl + X2 <^ X^ — 7- X2 + Xi 

2 

2 1 

X4 + X^ <^ Xq — Xl + X5 

1 



and the conjugacy constants ci = 1,C2 = 1,C3 = 2,C4 = 1,C5 = 2,cq = 1. 
This is easily identified as the enzyme network 

Si+E ^ SE ^ S2 + E 
S2 + F ^ PF ^ Si+F 

where an enzyme E facilitates the transfer of a substrate Si into another 
substract S2 and another enzyme F facilitates the transfer back. This net- 
work was considered extensively in p}] and |4j. In particular, it was shown 
in [4] that for all rate constant values, the network possesses within each 
invariant space a unique positive equilibrium concentration which is glob- 
ally asymptotically stable relative to that invariant space. It follows by the 



properties of linearly conjugate reaction networks that (29) inherits the same 
qualitative dynamics. 
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4 Structural Dynamical Equivalence 



In this section we extend the computation procedure given in Section [3] for 
dynamical equivalence to the case of networks M which are structurally fixed 
but have undetermined parameters. 



4.1 Dynamical Equivalence 

The MILP framework outlined so far requires that the rate constants be 
specified for the network M. Consequently, when we search for networks 
which are linearly conjugate to a given network M, we are really asking if 
there are networks which are linearly conjugate for a given choice of param- 
eter values. 

For networks where the dynamical behaviour is heavily dependent on 
the chosen rate constants, however, it is possible that certain behaviours are 
being overlooked by poor rate constant selection. There are networks, for 
instance, which are known to be linearly conjugate to weakly reversible net- 
works or complex balanced networks for certain values of the rate constants 
but not for others (see Examples 2 and 3 of [25j). In such cases, if the rate 
constants are not carefully chosen, the algorithm may overlook these net- 
works and we would not realize that the mechanism shares characteristics 
with these other networks. 

Therefore, we now change the problem setup by fixing only the structure 
of the initial network J\f but not the parameter values. We will show that 
this problem class can also be casted to the framework of MILP. This would 
remove the above mentioned limits of using a fully specified initial network 
model. 

The conditions for dynamical equivalence, keeping the entries of both 
Ak and A'^ general, are 



(DE) < 



Y-A'j^ = Y-Ak 



, m 



Y.^A'k^'^ = Y.^Ak]ij=0, 3 = 1,... 

1=1 i=l 

[A'k]ij > 0, [Ak]ij > 0, i,j = l,...,m, i^ j 
[A'k]ii<0,[Ak]ii<0, i 



(30) 



1, . . . ,m 



Note that in ( 30 ) both the off-diagonal entries of Aj. and are now decision 
variables. 



As in Section 2.3, we want to keep track of the structure of A',. The 
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conditions corresponding to ^ for the matrix A'j^ are 
f < [A[]ij - e6ij, ij = 1, 

(S2) 



< [A'i^]ij - e6ij, i,j = l,...,m, i j 
< -[A'klij +Uij6ij, i,j = l,...,m, i/j (31) 
6ij€ {0,1}, i,j = l,...,m,iy^j. 



As before, the binary variables 6ij keep track of whether a reaction is in the 
network J\f' or not and thus are capable of counting the number of reactions 
inAf'. 

We also, however, want to permit Aj. to have a variable rate constant 
values within a fixed network structure. In order to fix this network struc- 
ture, we introduce the binary variables 'yij G {0, 1}, i, j = 1, . . . ,m, i ^ j, 
and the logical equations 



7ij 



1 



[Ak]ij > e, i,j = l,...,m,i^j 



(32) 



for some < e ^ 1. In other words, the Ji/s keep track of whether the 
reaction Cj — )■ Ci is in the network M. The conditions required to allow the 
entries of A/, to vary independently within this pre-defined structure are 



(Ind) 



< 



k\ij ^lij^ 



1, 



< -[Ak]ij + Uij'~iij, i,j = l,...,m, j 



where 



7ij 



1, 
0, 



if {C„Ci)en 
otherwise. 



(33) 



(34) 



Further constraints can be implemented to search through subspaces of the 
parameter spaces for alternative realizations. For instance, if we suspect the 
reaction rate for the reaction Cjj — ?■ Cjj is slaved to that of Cj^ — >• Ci^ we can 
add [^fcliiii = [^fc]j2j2 to the procedure, etc. 

The conditions (30), ( |31[) a nd (33) can be combined with the structural 
conditions for reversibility (15) and weak reversibility (13 and 14) and the 



objective functions ^ and (10) to search over the parameter space of Af 
for sparse and dense alternative realizations A/"' which satisfy these further 
structural constraints. 



4.2 Complex Balanced Realizations 

It is also desirable to explore the parameter space of M for alternative re- 
alizations J\f' which are complex balanced. The linear constraints (22) and 
(18), however, cannot be used in the parameter-independent case since the 
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required equilibrium concentrations x* G M"q depend on the rate constants 
for N which are not specified. 

We might expect, however, since all weakly reversible networks are com- 
plex balanced for some choice of rate constants, that if a network N has 
a weakly reversible alternative realization N' for some other choice of rate 
constants then it also has a complex balanced alternative realization N" for 
some choice of rate constants. That is to say, in the parameter-independent 
optimization procedure weak reversibility is sufficient to demonstrate com- 
plex balancing. The main result of this subsection guarantees this (Theorem 



4.3D . 

First, however, we need the following results about weakly reversible 
networks. 

Theorem 4.1 (Theorem 3.1, [To]; Proposition 4.1, [16|). Let Ak he a kinetics 
matrix and let Aj, i = 1, . . . denote the support of the i^^ linkage class. 
Then the reaction graph corresponding to Ak is weakly reversible if and only 
if there is a basis of ker{Ak), {b^^^ . . . , b^^)}, such that, for i = !,...,£, 



b« 



> 0, j G A. 

= 0, j A, 




Theorem 4.2 (Theorem 1, flS^). Under the assumption of mass-action 
kinetics, weakly reversible chemical reaction networks possess at least one 
positive equilibrium concentration within each positive invariant space of the 
system. 



An immediate consequence of Theorem 4.1 is that a network is weakly 
reversible if and only if there is a vector b G M"o in the kernel of A^. We 
will exploit this fact in the next result. 

Theorem 4.3. Suppose there is a choice of rate constants such that the 
network M is dynamically equivalent to N' and N' is weakly reversible. 
Then there exists a choice of rate constants such that the network N is 
dynamically equivalent to M" where M" is complex balanced. Furthermore 
N" can be selected to have the same structure as N' . 

Proof. Let Ak be the kinetics matrix associated with N and be the 
kinetics matrix associated with AA', and suppose that N' is weakly reversible. 
Let b G M"o denote the positive vector in \.ei{Ak) guaranteed to exist by 
Theorem 4.1 and x* G M"q be any positive equilibrium concentration of ([T 



guaranteed to exist by Theorem 4.2 
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We now define a new network N" with the associated kinetics matrix 

(35) 



A'L = A'l, ■ diag j ^ 



where we define vector division to be componentwise division, i.e. x/y = 
[xi/yi,X2/y2, ■ ■ ■ ,Xn/yn] for X, y G M". Notice that all the terms in the 
definition of A'^ can be determined under the assumption that M' is weakly 
reversible. We have that 

<-M/(x*) = ^',.b = 

since b G ker(A^) so that M" is complex balanced at x* G IK>o- Furthermore, 
we have that 



Y ■A'l = Y -A'j,- diag — — }=Y-Ak- diag 



b 



^(x*)J ''[^(x*) 

It is clear that A^- diag{b/^'(x*)} has the same structure as Ak so that this 
corresponds to different choice of rate constants for the network M. The 
result follows. □ 

Note that, while we do not have linear constraints capable of determining 
complex balanced networks explicitly, we can always construct a complex 
balanced network from a weakly reversible or reversible network output by 
the optimization procedure. (Although this may not hold if further restric- 
tions on the parameter space of J\f have been imposed.) 

It is worth noting that the corresponding result for reversible and de- 
tailed balanced networks does not follow in the same manner as the proof of 



Theorem 4.3 This is because, for reversible networks with cycles, it is known 



that the detailed balancing condition entails further conditions on the rate 



constants above and beyond complex balancing 14, 18 . We can still, hoW' 



ever, construct a complex balanced realization from an arbitrary reversible 



network according to Theorem 4.3 Since detailed balancing implies no fur- 
ther dynamical information above and beyond complex balancing, for all 
practical purposes this is as far as we need to go. 

It is also worth noting that this algorithm cannot determine networks 
which are linear conjugate to a given structurally- fixed network; it can only 
find dynamically equivalent networks. This is clear since linear conjugacy 
requires that 

Y ■Ak = T-Y -Ab (36) 
where T = diagjc} and c G K™q. Regardless of whether we consider the 



transformation T on the left-hand-side or right-hand-side of ( 36 ) , it produces 
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a non-linear condition and therefore cannot currently be placed within the 
existing MILP framework. 



4.3 Examples 

In this section, we introduce a few examples which illustrate how the al- 
gorithm for producing weakly reversible, reversible, and complex balanced 
realizations of a structurally- fixed but parameter- variable network works. 



Example 3: Consider the reaction network N given by 



2Xi + X2 



3X1 



3X2 



Xi + 2X2 



where a > 0. Despite the reversible step in the central reaction, N is neither 
fully nor weakly reversible and therefore the dynamics do not fall within the 
scope of the theory of such networks. 

We want to check whether there are weakly or fully reversible networks 
which are dynamically equivalent to N for some value of a. We set Ci = 
2X1 + X2, C2 = 3X1, C3 = 3X2, and Ca = Xi + 2X2. Searching for 
a sparse network Q in GLPK over the weakly reversible constraints (30), 
dsTj ), ( pS] ), (13) and (14), with the additional constraints [Aa;]21 = [^fc]34 and 
[^fc]23 = [^fe]32 = 1 and the bounds e = 1/20 and u-ij = 20 i,j = 1,. . .4, 
i ^ j, gives the alternative realization 



2X1 + X2 ^ 3X1 



M' : 



3/2 I 
3X2 



1/20 



i3/2 
Xi + 2X2 



and a = 1/20 for the original network M. The network Af' has the corre- 
sponding kinetics matrix 



Ai 



__\_ 

20 








_3 
2 



3 
2 



3 
2 



_3 
2 
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and the positive equilibrium concentration (x*,X2) = (1, 1). It can easily be 
checked that M' is not complex balanced at this equilibrium concentration. 
In order to construct a complex balanced network J\f" with the same 



structure as J\f' by (35 ) we need to determine a vector b G Mt.Q such that b £ 
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ker(^'^). It can be easily checked that the vector b= [3/2 1/20 1/20 3/2]^ 
works; however, using this choice of b produces a set of rate constants by 



Ak ■ diag 



which clearly violates the condition [^fc]23 = [^fc]32 = 1- This can be solved 
by choosing another multiple of b. In fact, we can see that for the vector 
b = [30 1 1 30]"^ the appropriate rate constant choices for J\f occur by 
choosing a = 3/2 and the corresponding complex balanced realization given 



by (35) is 



2Xi + X2 

3/2 I 

3X, ^ 



3/2 



3Xi 

4-3/2 
2X2 



(This corresponds to a scaling by 3/2 of the 'block' network given in |24| . 
From the analysis presented in that paper, we know the network is complex 
balanced only for this particular value of a.) 

We may also be interested in fully reversible alternative realizations of 



N. Replacing the constraints (13) and (14) in the above algorithm by (15) 



and running in GLPK gives the network 

1/20 

M' : 2X1 +X2 ^ 3X1, 
3 



1/20 

X1 + 2X2 ^ 3X2 

3 



and a = 1/20 for the network N . The corresponding kinetics matrix is 



__!_ 

20 




















j_ 
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which has the positive equilibrium concentration (x*,X2) = (1,1) which is 
neither complex nor detailed balanced. In order to find a complex balanced 
network with the same structure as M' we notice that the kernel of A'^ is 
given by the span of [3 1/20 0]^ and [0 1/20 3]^. In order to 
preserve the property [^fe]23 = [^A;]32 = 1 for the network M we need to 
choose b = [60 1 1 60]'^. This gives the value of q = 3 for the network N 
and the complex balanced realization 



M" : 2Xi+X2^3Xi, 
3 



Xi + 2X2 «^ 3X2. 

3 
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It can easily be checked that N" is also detailed balanced. 
Example 4: Consider the reaction network N given by 

N : 2X1 ^ 2X2 
and the reversible alternative realization N' given by 



1/2 1/2 1/8 

2X1 ^ X1+X2 T± 2X2 T± 2X1. 

1/2 3/4 1/4 



If we make the associations Ci = 2Xi, C2 = Xi + X2 and C3 
then the network Af' has the corresponding kinetics matrix 



2X0 



Ai 



lyQ we have A'f^ 



^(x*) / 



Choosing any equilibrium concentration x* G 
so that the network is not complex balanced. 

We wish to construct a complex balanced realization M" using the method- 
ology outlined in the proof of Theorem 4.3 We choose the equilibrium 
value ixl,X2) = (1,1) so that *(x*) = [1 
b = [1 5/4 1]"^. According to (35) we have 



1 1]-^ and the kernel vector 





' 1 





■ 









5 
4 















1 





It can be easily checked that the corresponding network M" is complex 
balanced and is dynamically equivalent to A^. (In general the rate constants 
of J\f may change; however, in this case we have A/.- diag{b/^'(x*)} = A^.) 
Despite being complex balanced and reversible, the realization M" is not 



detailed balanced according to (20). We might wonder if we can construct 
an alternative realization in a similar manner as (35). We can easily see, 
however, that applying the detailed balancing conditions to Ak- diagjc}, 
where c E M5,q is arbitrary, produces the unsatisfiable set of conditions 



Cl 2 



C2 



-X1X2, 



Cl 2 

4^1 



C3 2 
77^25 



C2 



-X1X2 



3C3 

4 



X2- 



In other words, we cannot always construct a detailed balanced network 
from an arbitrary reversible network M' in the same constructive manner 
used for complex balanced networks from weakly reversible networks. 
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5 Conclusions and Future Work 



In this paper, new computational methods were presented for finding net- 
works which are hnearly conjugate and dynamicahy equivalent to a given 
chemical reaction network endowed with mass-action kinetics. 

It was demonstrated that finding a dense or sparse reversible, detailed 
balanced and complex balanced network which is linearly conjugate to a 
given network can be framed as a MILP optimization problem. The case of 
determining conjugate networks which have the greatest and fewest number 
of complexes was also extended to the case of non-trivial linear conjugacies. 

It was shown that, similarly to the case of dynamical equivalence j37j, the 
graph structure of linearly conjugate dense networks containing the maxi- 
mal number of nonzero reaction rate coefficients is unique, and that the 
unweighted directed reaction graph of any linearly conjugate network is a 
proper subgraph of the unweighted directed reaction graph of the dense 
one if the set of complexes is given. Moreover, it was proved that arbitrary 
equilibrium points of the initial network can be used for the existence check- 
ing and computation of linearly conjugate complex balanced and detailed 
balanced networks. 

Additionally, the problem of dynamical equivalence was studied when 
only the structure of the initial network N is fixed, but its rate constants 
can take any positive value. It was shown that in this case, the computation 
of dense and sparse weakly reversible and reversible networks M' can also be 
formulated as a MILP problem. It was also shown that complex balanced 
networks M" with identical structure to N' can be constructed from these 
realizations and corresponded to an alternative choice of rate constants for 
N . This modification allows us to scan through a range of parameters 
as part of the procedure and therefore answer questions about a network 
based on its structure alone. The operation of the developed methods were 
illustrated on numerical examples. The achievements further extend the 
applicability range of many existing structure-dependent results of chemical 
reaction network theory. 

Further areas of research and open questions include: 

1. The procedure introduced to determine structurally dynamically equiv- 
alence networks is as of yet unable to search through non-trivial lin- 
early conjugate networks in a linear manner. As such, many networks 
with potentially insightful information about the dynamics of a given 
network are being overlooked. Incorporating non-trivial linear conju- 
gacy into a manageable optimization framework is therefore of primary 
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interest. 



2. The results obtained so far depend on the reaction networks under 
consideration having mass-action kinetics ([T]). Conjugacy and compu- 
tational results for other widely-used kinetic schemes (e.g. Michaelis- 
Menten kinetics |29|, Hill kinetics (2ll) would greatly expand the scope 
of applicability of these methods. 

3. There are many classes of networks with known behaviour lying out- 
side the scope of weakly reversible network theory (4][8j[9 . Determining 



constraints which could restrict our search to within these classes of 
networks would broaden the scope of dynamical behaviours (e.g. peri- 
odicity, oscillatory behaviour, multistability, etc.) we could guarantee 
through this computational procedure. 
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